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Abstract 


The Yee finite difference time domain (FDTD) algorithm is widely used in computational electromagnetics because of its 
simplicity, low computational costs and divergence free nature. The standard method uses a pair of staggered orthogonal 
cartesian meshes. However, accuracy losses result when it is used for modelling electromagnetic interactions with objects 
of arbitrary shape, because of the staircased representation of curved interfaces. For the solution of such problems, we 
generalise the approach and adopt an unstructured mesh FDTD method. This co-volume method is based upon the use of a 
Delaunay primal mesh and its high quality Voronoi dual. Computational efficiency is improved by employing a hybrid primal 
mesh, consisting of tetrahedral elements in the vicinity of curved interfaces and hexahedral elements elsewhere. Difficulties 
associated with ensuring the necessary quality of the generated meshes will be discussed. The power of the proposed solu- 
tion approach is demonstrated by considering a range of scattering and/or transmission problems involving perfect electric 
conductors and isotropic lossy, anisotropic lossy and isotropic frequency dependent chiral materials. 


1 Introduction 


In this paper, we consider the simulation of problems in 
electromagnetics using unstructured co-volume staggered 
meshes and a generalisation of the Yee finite difference 
time domain (FDTD) method [1]. In its standard form, the 
FDTD algorithm is implemented on structured staggered 
spatial meshes. The use of staggered meshes enables the 
interleaving of the unknown nodal values of the electric and 
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magnetic field components, resulting in an explicit algorithm 
that is second order accurate in both space and time. The 
algorithm has the additional advantages of simplicity, low 
computational cost and, as no linear algebra is required, 
there is no inherent limit to the size of the simulation. In 
addition, the scheme preserves the energy and amplitude 
of waves, is divergence free and can readily be parallelised. 
To retain these benefits, for simulations involving geom- 
etries of complex shape, non-uniform and unstructured mesh 
FDTD implementations have been suggested, such as the 
generalised Yee algorithm [2] and the Yee-like algorithm 
[3]. Unfortunately, these methods do not retain the effi- 
ciency of the original scheme [4]. An alternative approach 
to generalising the FDTD algorithm to unstructured meshes 
is to employ a primal unstructured Delaunay mesh and its 
orthogonal Voronoi dual graph [5]. With this method, hybrid 
meshes can be naturally incorporated, without requiring the 
inter-mesh interpolation, and consequent non-physical dif- 
fraction effects mesh interfaces, that plagues many hybrid 
approaches. The proposed implementation is edge and face 
based, which means that it is not restricted to specific mesh 
element types and can handle any form of polyhedron. 
This paper is structured as follows: In Sect. 2, we outline 
the problem of interest and employ a scattered field for- 
mulation for Maxwells equation in integral form for prob- 
lems involving an isotropic lossy dielectric material and/ 
or a perfect electric conductor (PEC). Section 3 introduces 


val Springer 


182 


the FDTD method and illustrates the discretisation of the 
Maxwell equations on both a structured and an unstructured 
mesh. Section 4 is devoted to the mesh generation process, 
which is one of the most crucial parts for the successful 
implementation of the proposed scheme. The treatment of 
the different kinds of boundary conditions is discussed in 
Sect. 5. Section 6 describes the method adopted for the cal- 
culation of the distribution of the radar cross section (RCS). 
In Sect. 7, we validate our algorithm by comparing the 
numerical and analytical solutions for the RCS distribution 
for problems involving, dielectric lossy and coated spheres 
of different electrical lengths. The calculation of the RCS 
distribution of a more complex PEC aircraft geometry and 
the simulation of transmission of an electromagnetic pulse 
through a radome are presented in Sect. 8. In Sects. 9-11, 
we move on to consider anisotropic lossy materials, where 
the material parameters, electric permittivity ¢, the mag- 
netic permeability and the electric and magnetic conduc- 
tivites o,0,,, become second order tensors. This allows us 
to consider the simulation of problems involving crystals, or 
composite materials, where the orientation of fibres plays an 
important role. In Sects. 12-13, we consider a further exten- 
sion that allows for the modelling of bi-isotropic materials. 
Materials of this type are frequency dependent and have an 
additional parameter, the chirality, which couples the electric 
to the magnetic field, making the simulation even more chal- 
lenging. Brief conclusions are presented in Sect. 14. 


2 Problem Formulation 


We consider problems involving the simulation of the inter- 
action of an incident electromagnetic wave, generated by a 
source located in the far field, and a general medium sur- 
rounded by free space. We assume, initially, that the medium 
might represent a dielectric or a closed region of very high 
electrical conductivity, as illustrated in Fig. 1. A closed 
region of very high electrical conductivity is generally 
approximated as a PEC. With this approximation, electro- 
magnetic fields do not penetrate the conducting surface and 
the outer boundary of the PEC forms the inner boundary of 
the problem domain. For such simulations, it is convenient 
to split the total electric and magnetic fields into incident and 
scattered components, according to 


Ero: = Eine + E 


mc scat 


H = Hine + Hycat (1) 


tot 


with the subscripts tot, inc and scat referring to the total, 
incident and scattered fields respectively. To enable a 
numerical solution of this problem, the governing Maxwell 
equations are expressed in terms of the scattered field com- 
ponents [6] and written in the integral form provided by the 
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Dielectric or Perfect Conductor (PEC) 


Fig. 1 Schematic representation of the electromagnetic scattering 
problem 


Laws of Ampére and Faraday [5]. For a three dimensional 
lossy dielectric medium, the resulting equations become 


IE ine 
- f (e-) 3 OA (2) 


ot 
~ | aH “dA — [oF -dA 
A A 
(3) 


Here, 0A denotes the closed curve bounding a surface A, dA 
is an element of area of the surface A, directed normal to the 
surface, and dl is an element of contour length in the direc- 
tion of the tangent to the curve 0A. In addition, € denotes 
the electric permittivity of the medium, y its magnetic per- 
meability, o its electric conductivity and o,, its magnetic 
conductivity. 

We will often consider plane incident waves, of free space 
wavelength A. In this case, if €, and yz, denote the electric 
permittivity and magnetic permeability of free space respec- 
tively, we can write 

kxE,,, 


Einc = Eo cos(@t —k-r) H,,. = (4) 
q 


where k is the unit wave vector, r is the position vector, 


@ = 2z/A,is the angular frequency and 4 = 1/y,/é, is the 
impedance of free space. 
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3 The FDTD Method 
3.1 The Standard Algorithm 


The standard FDTD scheme [1, 7] is a popular solution 
algorithm in computational electromagnetics. In its origi- 
nal form, the algorithm is implemented on two mutually 
orthogonal staggered structured cartesian meshes in a 
straightforward manner. 

For example, on a general (x, y, z) cartesian mesh, with 
mesh spacing Ax, Ay, Az in the three coordinate directions, 
the location of the unknown electric and magnetic field 
components is illustrated in Fig. 2. These components are 
advanced in time using a staggered leapfrog scheme. In 
the resulting algorithm, with the superscript n denoting 
an evaluation at time nAf, where Af is the time step, and 
subscripts (/, J, K) denoting an evaluation at the point 
([Ax, JAy, K Az), typical discretized forms of Eqs. (2) and 
(3) in the scattered field formulation [6] are 
(Qe + CADE eevee 

= (Ze - CAN) E atx (l4! /2.J,K) 


n+1/2 
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and 


Fig. 2 Location of the unknowns for the cartesian FDTD algorithm 
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Similar equations can be written directly for the other field 
components. For a stable implementation of this explicit 
scheme, the magnitude of the time step employed has to be 
restricted and, for a uniform cartesian mesh, this stability 
criterion may be written as cV3At < 1, where c = (ey)~'/” is 
the speed of wave propagation through the medium and / is 
the edge length in the mesh [7]. It should be noted that, since 
the incident fields are defined analytically in Eq. (4), the time 
derivatives of E;,,,. and H;,,, can be evaluated directly for 
use in Eqs. (5) and (6). In this form, this explicit algorithm is 
second order accurate, in both space and time, on a uniform 
mesh. The computational implementation of this FDTD 
method, characterized by its low storage requirements, of 
around 75 Mbytes per million nodes, and its low operation 
count, provides a highly efficient simulation process. 


3.2 Unstructured Mesh Algorithm 


Although the cartesian mesh algorithm is very efficient, 
unstructured meshes are more suitable for use in the solu- 
tion of problems involving regions of complex shape. For 
the implementation of a corresponding unstructured mesh 
FDTD procedure, a suitable pair of mutually orthogonal 
meshes has to be identified. An obvious choice, given the 
widespread availability of automatic unstructured mesh 
generators, is to select a primal Delaunay tetrahedral mesh 
together with its Voronoi dual. These meshes are mutually 
orthogonal, in that each Voronoi face is a perpendicular 
bisector of the corresponding Delaunay edge, while each 
Delaunay face is perpendicular to the corresponding Voronoi 
edge. It is assumed that the Delaunay mesh employs NP and 
the Voronoi mesh NY edges. Each Delaunay edge has an 
associated closed loop of Voronoi edges and each Voronoi 
edge is surrounded by a closed loop of Delaunay edges. The 
unknowns in the unstructured mesh implementation of the 
FDTD scheme are located at the midpoints of these edges. 
The unknown at the centre of the ith Delauany edge corre- 
sponds to the projection, E,.,, ; of the scattered electric field 
onto the direction of the edge. The unknown at the centre of 


the jth Voronoi edge corresponds to the projection, Hy.4,;, of 
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the scattered magnetic field onto the direction of the edge. 
Discretization of the Ampére’ and Faraday’ Laws by apply- 
ing the central difference approximations 


Et —k" AY 
[ec dA & € (Excot ~ Evcar) AY 


5 Esca a 7) 
A 

MP 

ae n 1D 
¢ E, .,dl © 2, Ea (8) 
aA 7 
ntl n 4 
n+1/2 o( Ett + a )A 

/ OK, .dA oa eke AY = 2 (9) 
A 


wee? = a 


r scat scat 

—H,,,,dA & Mm 
/ Hae scat J\t : 
A 

MY 
os n+1/2 WV 

¢ Ayal x Dy As cat jy i oa 
0A ~ 


‘scat scat 


o( Hy” +H )A? 
/ On seqdh x OH" AP == 
A 


2 
(12) 
leads to the equations 
Qe+o ADE! =(Qe-0 ADE”, 
OE... n+1/2 
-2A te - &)) >= 
: (13) 
A MY 
2A tl] oO 41/2 Vv V pntl/2 
a AY 2 1: ee a oA; E inci 
and 
Qu to, Aditi? =(Qu-o, AdH 
OH ine j e 
—2 A tu— uy) a 
D 
2 A t n 1D 
+ D > » Pots a 
ij k=1 : 
_ Dyn 
Om; a (14) 


Here Ip represents the length of the ith Delaunay edge and 
a corresponds to the area of the Voronoi face spanned by 
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the Voronoi edges surrounding Delaunay edge i. Similarly, 
r represents the length of the jth Voronoi edge and a? cor- 
responds to the area of the Delaunay face spanned by the 
Delaunay edges surrounding Voronoi edge j. The numbers 
Jig k = 1,...,MY/ refer to the MY edges of the Voronoi face 
corresponding to the ith Delaunay edge, as illustrated in 
Fig. 3a. Similarly, the numbers lies o_o nee ,»M? refer to the 
M? edges of the Delaunay face corresponding to the jth 
Voronoi edge, as illustrated in Fig. 3b. 

With these staggered equations, the magnetic field is 
advanced over the dual mesh at the half time step, using 
Eq. (14), and the electric field is advanced over the primal 
mesh at the full time step, using Eq. (13). It should be noted 
that, with the appropriate interpretation, the discrete approxi- 
mations of Eqs. (13) and (14) reduce to the basic FDTD algo- 
rithm, when they are employed on a pair of staggered regular 
orthogonal cartesian meshes. For post-processing purposes, 
a local least squares problem has to be solved to determine 
all components of E and H at the vertices of the primal and 
dual meshes respectively. For a stable implementation of this 
explicit unstructured mesh scheme, practical experience has 
shown that a relationship of the form 


cAt < S, min{/;1?} (15) 
generally results in stability, where S; is a safety factor that 
depends on the mesh and is normally given a value in the 
range 0.8-2.0. 

4 Mesh Generation 


4.1 Structured Mesh 


Structured meshes of hexahedra are readily generated and 
the meshes employed for wave propagation problems are 


Fig.3 a The ith Delaunay edge, connecting Delaunay vertices p, and 
P», and the corresponding Voronoi face, formed by the Voronoi edges 
Jit» +++ >Jig; b the jth Voronoi edge, connecting Voronoi vertices e, 
and e,, and the corresponding Delaunay face, formed by the Delaunay 
edges i; ,i;2, 113 
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generally constructed to consist of cubes of uniform side 
length, 6. For problems involving waves of characteristic 
wavelength, Ap, a value of 6 in the range A) /30-A,/10 is 
typically employed for practical applications of the FDTD 
scheme on regular cartesian grids [7]. 


4.2 Unstructured Mesh and the Ideal Unstructured 
Mesh 


The unstructured mesh discretisation employed in Eqs. (13) 
and (14) will be second order accurate provided that each 
Voronoi face is a perpendicular bisector of the correspond- 
ing Delaunay edge. This criterion is satisfied automatically 
for meshes generated by the Delaunay method, but the need 
to recover a given boundary triangulation may result in the 
loss of this property. In addition, each Delaunay face must 
be a perpendicular bisector of the corresponding Voronoi 
edge and the intersection point, between the edge and the 
face, must be at the centroid of the face. This criterion can 
only be achieved for a perfect tetrahedron, which is a tetrahe- 
dron for which all faces are equilateral triangles. In this case, 
the centre of the circumsphere lies inside the tetrahedron 
and coincides with its centroid. Failure to ensure that this 
criterion is satisfied means that the corresponding contour 
integral in Eq. (13) is approximated over a surface that does 
not intersect the corresponding Voronoi edge, as illustrated 
in Fig. 4. In this case, the accuracy of the approximation of 
the integral cannot be guaranteed. 

It follows that certain requirements should be imposed 
on a Delaunay primal mesh and its Voronoi dual mesh, to 
ensure the successful implementation of the unstructured 
mesh FDTD solution method. These requirements are: 


1. all Voronoi and Delaunay edge lengths should be 
bounded from below; 

2. the centre of the circumsphere for each Delaunay tetra- 
hedron should lie within the tetrahedron; 

3. all Delaunay and Voronoi edge lengths should be 
bounded from above by a value that is not significantly 
greater than the specified value of 6; 

4. any deviation in the location of the midpoint of a Voro- 
noi edge from the actual point of intersection with the 
corresponding Delaunay face should be minimised; 


Fig. 4 Illustration of the case 
where a Voronoi edge and the 
corresponding Delaunay face do 
not intersect 
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5. any deviation in the location of the circumcentre of a 
tetrahedron from its centroid should be minimised. 


In this list of requirements, the first two entries are the most 
important in securing a practical, stable implementation of 
the unstructured FDTD algorithm. For complex three dimen- 
sional geometries, it is not always possible to satisfy exactly 
the last two requirements, so that these criteria have to be 
relaxed in practice. The third requirement guarantees a suf- 
ficiant resolution for the propagation of the electromagnetic 
wave. 

In two dimensions, the ideal unstructured mesh is the 
space filling mesh of equilateral triangles [8]. The number 
of nodes connected to each node is 6, the Voronoi edge 
length l” = 2 /1/3 ~ 0.561 and all the vertex angles are 
60°. In three dimensions, a tetrahedralmesh made of ele- 
ments with faces in the form of equilateral triangles is not 
space filling. It has been shown [9-11] that an ideal space 
filling mesh consists of equal non-perfect tetrahedra, with 
each face in the form of an isosceles triangle, with one side 
of length Lo and two shorter sides of length 


Pt 3/2 ae ng’ SiX Such tetrahedra form a parallelepiped 
tiling the space, as illustrated in Fig. 5. This configuration 
can be shown to maximise the minimum Voronoi edge for a 
fixed element size. All Voronoi edges have the same length 
IY = 0.38 d6whered = (1?) = (3 fe + 4[0 /T = 0.92 Feng’ 
The number of nodes connected to each node is 14. All ver- 
tices have an acute angle of 71.5° and, hence, for each ele- 
ment, the centre of the circumsphere is located inside the 
element. However, certain dihedral angles are equal to 90°, 
which is larger than the value 70.5° for the dihedral angle of 
the perfect tetrahedron. The elements in this three dimen- 
sional ideal mesh satisfy the requirements set out above, but 


do not, of course, fit general boundaries. It can be readily 


Fig.5 Detail of a mesh of ideal tetrahedral elements, showing the 
surface Delaunay faces and the internal Voronoi cells 
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demonstrated that the solenoidal preserving property of the 
original FDTD algorithm is also valid for the discretisation 
of Eqs. (13) and (14) on these ideal meshes. 


4.3 Unstructured Mesh Generation 


Traditional automatic unstructured mesh generation meth- 
ods, such as the advancing front technique [12] and the 
Delaunay triangulation [13], or their combination [14], 
are not designed to guarantee the creation of an unstruc- 
tured mesh meeting the FDTD requirements set out above, 
even in the two dimensional case. These methods gener- 
ate meshes in which the element edge length is accept- 
able, but the corresponding Voronoi diagram is often 
highly irregular and can include some very short Voronoi 
edges. This means that these methods cannot guarantee 
the regularity of the edge lengths of the dual mesh and 
the absence of bad elements. Although methods based 
on swapping, reconnection and smoothing [15] can be 
used to improve mesh quality, experience has shown that 
boundary constraints usually mean that their use does 
not result in a suitable unstructured mesh for the FDTD 
method, even in two dimensions [8]. 

An alternative approach is to generate meshes using 
an iterative constrained centroidal Voronoi tessellation 
(CVT) [16], followed by mesh quality optimisation. The 
iteration process adopted is Lloyd’s algorithm [17] which, 
given an initial mesh, relocates the nodes, to the mass 
centroids of the corresponding Voronoi cells, and then 
creates a new Voronoi tessellation of the relocated nodes. 
The process is repeated until all nodes are close enough to 
the corresponding centroids. The initial tetrahedral mesh 
is constructed using the vertices of the ideal mesh as the 
point distribution. In addition to relocating the nodes, the 
CVT scheme changes the mesh topology and, although 
the quality of final mesh is much higher than the quality 
of the initial mesh, typically around 10% of the elements 
in the final mesh may still be found to violate the mesh 
quality measures that are required in the current context. 


4.4 Unstructured Mesh Optimisation 


The unstructured mesh is optimised in an attempt to ensure 
that both the primal and the dual mesh are of the highest 
possible quality. The requirement that a dual edge must be a 
bisector of the corresponding Delaunay edge is relaxed. At the 
same time, the corresponding dual mesh vertex is moved to a 
point which still ensures orthogonality between the two grids 
and which lies inside the corresponding primal element. The 
procedure adopted is similar to that found in the concept of a 
power diagram [18], where the perpendicular bisector of an 
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edge is moved parallel to itself by an amount that is propor- 
tional to a weight assigned to the edge nodes. 

To illustrate how this optimisation process can be 
achieved, consider, in two dimensions, the case when the 
Voronoi vertex, located at the centre, O, of the circumcircle 
enclosing a triangular element ABC, lies outside the element. 
The point O is located at the intersection of the perpen- 
dicular bisectors of the three edges of the triangle. For edge 
AB, the perpendicular bisector is identical to the common 
cord of the two circles, each having the same radius as the 
circumcircle, centred at the points A and B, as illustrated in 
Fig. 6a. When the radii of these two circles are changed, the 
new common cord remains perpendicular to AB but loses 
its property of being a bisector. This is illustrated in Fig. 6b. 
In this case, with careful choice of the radii, the common 
cords of the circles centred at the nodes of each edge still 
intersect at a common point O. The vertices and edges of 
the dual mesh are modified accordingly, with the dual mesh 
vertex for the element being relocated to this common point. 
This same process can be applied to a tetrahedral element in 
three dimensions. 

These ideas are employed as the basis of a mesh optimi- 
sation process that is used to locate the position XQ of the 
dual mesh node, O, for each element, as close as possible 
to the element centroid. This will ensure that each element 
contains its dual node and guarantees the essential criterion 
of mesh orthogonality. For a tetrahedron element, ABCD, 
the dual node is located at the intersection point of the com- 
mon cords of the spheres centred at the vertices of each 
element, with the sphere at A having a radius of R,. If the 
radius employed at each vertex of a tetrahedron is the same, 
then the dual node will coincide with the circumcentre. Sup- 
pose that Xy = (x°, y?, z°) and that the vertex coordinates 
are X, = (x", y¥, z®) for E = A, B, C, D. The optimisation is 
achieved, during the CVT process, by changing the radius 
associated with vertex A by an amount 6R,, where the new 
dual node location for each tetrahedron is computed by solv- 
ing the system 


(i? — iA) — 6R,? = G? — i? — 6R? for E=B, Cand D 


(16) 


where the summation convention is employed for i = x, y, z. 
4.5 AHybrid Mesh Approach 


The unknowns in the FDTD algorithms are associated 
with the edges of the primal and dual meshes. It fol- 
lows that, for a fixed number of points per wave length, 
an unstructured mesh made of ideal tetrahedral elements 
will over-resolve the wave, by about 30%, compared to a 
structured cartesian mesh. To obtain a resolution simi- 
lar to that achieved with a structured cartesian mesh, the 
number of points per wavelength must be reduced when 
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(b) 


Fig.6 Mesh optimisation: a a two dimensional element ABC, with the centre, O, of the circumcircle lying outside the element, b a two dimen- 
sional element ABC, with the dual vertex moved to the point O inside the element 


generating the ideal tetrahedral mesh. Unfortunately, such 
a reduction can result in the under-resolution of bound- 
ary geometries, particularly in regions of high curvature. 
This suggests that a hybrid mesh approach could result in 
a computationally effective compromise. This has been 
demonstrated for two dimensional electromagnetic scat- 
tering simulations, which have been performed accurately 
and efficiently using such a hybrid approach, with a mesh 
of triangles used near complex boundaries and a cartesian 
mesh used elsewhere [19]. Initial attempts at extending 
these ideas to three dimensions, for problems involving 
simple geometries, have also been reported [11, 20]. 
With such an approach, the resolution employed is 
selected to be that appropriate for the operation of the 
FDTD algorithm on a regular cartesian mesh. This will 
result in an unstructured mesh that over-resolves the solu- 
tion in the free space region immediately adjacent to the 
scatterer. However, as it can be estimated that, typically, 
only around 10% of the total number of nodes will lie in 
the unstructured mesh region in three dimensions, follow- 
ing this approach will only result in a small increase in the 
total number of degrees of freedom, which is regarded as 
being acceptable given the improved boundary resolution. 


4.6 Hybrid Mesh Generation 


Consider the problem of discretising a general domain, using 
tetrahedral cells in the vicinity of the scatterer and regular 
cartesian hexahedral cells elsewhere. To provide a consist- 
ent hybrid mesh, these two groups of cells will be connected 
through a set of interface polygonal cells. The process of 
generating a hybrid mesh for the computational domain sur- 
rounding a general scattering obstacle is accomplished in 


four stages. In the first stage, an unstructured triangulation 
of the surface of the scatterer is produced [21] and the tri- 
angulation is then placed inside a hexahedral box. The sur- 
face of the box forms the outer surface of the computational 
domain. The region inside the box is discretised using a reg- 
ular cartesian mesh of cubes, of a prescribed edge length 6, 
as illustrated in Fig. 7a. Cubes within a prescribed distance 
of the scatterer, or lying internal to the scatterer, are removed 
in a second stage, to create a staircase shaped surface that 
completely encloses the scatterer, as shown in Fig. 7b. In 
the third stage, a point distribution is specified to completely 
cover the unmeshed region, as in Fig. 7c, with those points 
from the distribution that lie either inside the scatterer or 
outside the staircase surface being removed, as shown in 
Fig. 7d. The fourth stage consists of using the points that 
remain to generate an unstructured tetrahedral mesh in the 
region between the surface triangulation of the scatterer and 
the staircase surface. However, this mesh generation process 
cannot start directly, as the exposed faces of the staircase 
surface are all squares. The traditional approach to handling 
the problem of fitting a hexahedral mesh to a tetrahedral 
mesh is either to place a pyramidal element on each exposed 
square face, leading naturally to a consistent mesh, or, at 
the expense of introducing a hanging edge, to divide each 
exposed square face into two triangles. The use of pyramids 
is not adopted here, because the geometrical constraints of 
the staircasing imply that the centre of the circumsphere 
of each pyramid would lie outside the pyramid, which is 
something that should be avoided. Instead, the hanging edge 
approach is employed, but is modified to produce a consist- 
ent mesh of arbitrary shaped polyhedral cells that are appro- 
priate for use with the co-volume method. In this approach, 
each cube that forms part of the staircase surface is divided 
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Fig. 7 Illustration of the stages 


involved in the generation of 


an initial hybrid mesh: a the 


regular cartesian mesh covering 


the interior of the computa- 


tional domain; b creation of the 


staircase surface at a prescribed 


distance from, and completely 


enclosing, the scatterer; c 


the introduction of the point 


distribution; d the final point 


distribution, retaining only 


those points lying in the region 
between the scatterer and the 


staircase surface 


into six tetrahedra. For a cube with a single exposed face, a 
tetrahedron that contains one of the staircase surface trian- 
gles is removed from the mesh, as shown in Fig. 8a and b. A 
similar procedure is employed for cubes with either two or 
three exposed faces, producing the configurations shown in 
Fig. 8c and d. This process will replace the original staircase 
surface by a set of triangular faces that are appropriate as 
the starting point for the generation of the tetrahedral mesh. 
Each polyhedral interface cell is stored as a set of tetrahedra. 
The tetrahedra in each set have the same circumsphere and 
they are merged again during the solution process. 


(a) (b) 


4.7 Element Merging 


The primal Delaunay mesh will be degenerate if two, or 
more, adjacent tetrahedral elements share the same circum- 
sphere. In this case, certain Voronoi edges will have zero 
length, certain Voronoi faces will have zero area and the 
discrete approximations of Eqs. (13) and (14) become inva- 
lid. However, the approximations will remain valid if, at 
the outset, each set of such adjacent degenerate tetrahedra 
is merged, to create a polyhedron of the appropriate shape, 
and the unknowns and the lists of Delaunay and Voronoi 
edges are modified accordingly. The primal mesh is then 
regarded as a hybrid of tetrahedra and polyhedra and the 
dual mesh is updated appropriately. To avoid creating a new 


(c) (d) 


Fig.8 Treatment of the cubes forming the staircase surface to enable a consistent hybrid mesh to be generated: a typical cube; b cube with one 
face exposed; ¢ cube with two exposed faces; d cube with three faces exposed 
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data format to cater for the shapes that are created in this 
way, a polyhedron is stored as a set of tetrahedra, each with 
the same circumsphere. These tetrahedra are merged again 
during the solution process. 

Voronoi edges having a ’small’ length may still be present 
in the mesh and this will have a detrimental effect on the 
efficiency of the explicit time stepping process. To alleviate 
this problem, certain ’small’ Voronoi edges are also col- 
lapsed and the corresponding primal elements are merged 
to form a polyhedron. For the examples that are presented 
here, this merging is performed for all Voronoi edges with a 
length that is less than 0.016. Note that, if two merged poly- 
hedra are adjacent, it is important to ensure that the triangles 
forming the common face are co-planar. Violation of this 
criterion implies that the original merging process should 
not be allowed. Here, two triangles are taken to be co-planar 
provided their dihedral angle in less than 1°. 


5 Treatment of Boundary Conditions 
5.1 Material Interfaces 


At interface boundaries, the material parameters change 
over the area of integration. In this case, the values of €, py, 
o and o,, are replaced by appropriately averaged values €,,,, 
Hay» Sg, and o,,,,. For the standard FDTD algorithm, the 
three most common averaging techniques, at an interface 
between materials with generic parameters a, and da), are 
the arithmetic mean, the harmonic mean and the geometric 
mean. It has been reported [22] that the convergence rates 
obtained with the geometric and harmonic means dete- 
riorate as the mesh size decreases, whereas second order 
convergence is maintained with the arithmetic mean. For 


this reason, the arithmetic mean average is adopted here 


Fig.9 Averaging of material 

properties at an interface. a 

Delaunay edges in material | 

in blue, in material 2 in green, 

at the interface in red; b the 1,01 
Voronoi edge crossing a triangle 

at this interface in black. (Color 

figure online) 
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and modified for use with unstructured meshes, where the 
cell size is not necessarily uniform. 

To illustrate how an arithmetic mean is employed in 
the construction of the material parameters, consider the 
situation at an interface, between material 1 and material 
2, illustrated in Fig. 9. At the interface, material property a 
on Delaunay edge i is assigned the averaged value a 
computed as 


avDel,i? 


2M vor 
a Wr cell(k) 
CT alae a (17) 


In this equation, the term a,,,74, denotes the material param- 
eter assigned to cell k, surrounding Delaunay edge i. Cell k, 
with volume w,, is defined as the region spanned by the end 
points of Delaunay edge i, the point of intersection of the 
Voronoi edge with the Delaunay face and the position of the 
circumcentre of cell k. For example, in Fig. 9a, w, would be 
the volume spanned in material 1 by the points P1, P5, P9, 
P8, and a.gy(4) = € OF Aeeiy1) = 0. Similarly, w, would cor- 
respond to the volume spanned in material 2 by the points 
P1, P5, P9, P10, so that 4,42) = €2 OF Acei(2) = F2- 

Similarly, material property b, on a Voronoi edge, j, is 
assigned the average value b,,,,,; and is computed using 
the weighted average formulae 


2 
_ a BoD ceil) 


avVor,j —_ 2 
24 &¢ 


In this equation, the term b_,,,(¢) denotes the material param- 
eter assigned to cell 7. The length of the Voronoi edge inside 
cell 7 is g, and this corresponds to the distance between the 
point of intersection of the Voronoi edge j with the Delau- 
nay face and the circumcentre of the cell. For example, in 


(18) 


Hy, Om1 


(b) 
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Fig. 9b, the distance between the points P1 and P3 is g,, and 
the coefficient b,ey1) = Hy OF Deen) = Om, 

A PEC boundary may also be simulated in this fashion, 
treating it as a resistive sheet. In this case, the Delaunay 
and Voronoi edges forming the interface are assigned a 
very high conductivity value in the order of 10°S/m. 

It can be readily demonstrated that this averaging process 
on an unstructured mesh corresponds to the arithmetic mean 
on a structured mesh [23]. 


5.2 Far Field Boundary 


In the far field, the scattered fields should consist of outgo- 
ing waves only. This radiation condition is imposed by con- 
sidering the outermost layers of the structured hexahedral 
mesh to be in the form of an artificial absorbing perfectly 
matched layer (PML) [24]. The PML region is defined to be 
of a constant thickness and a standard PML formulation is 
adopted [7, 25]. 


6 Calculation of the RCS 


The RCS provides a measure of the power scattered from a 
target obstacle towards a receiver. In numerical simulation, 
we evaluate the bi-static RCS, which is a function of aspect 
angle and bi-static angle. The scattered wave is decomposed 
into two components. The first component has the same 
polarization as the transmitted signal and is referred to as 
co-polarized, while the second component has an orthogo- 
nal polarization state compared to the transmitted signal 
and is referred to as cross-polarized. For numerical simula- 
tions, the RCS distribution is of practical interest because an 
analytical solution is available for certain simple cases. To 
compute this distribution, a near to far field transformation 
procedure [26] is employed to obtain the required far field 
information from the computed near field data. To achieve 
this, a closed collection surface, S, completely enclosing the 
scatterer, is constructed and, when steady state conditions 
have been achieved, a further cycle is computed. During this 
cycle, the harmonic solution produced by the time domain 
solver is used to calculate the phasors, which enable fic- 
titious electric and magnetic currents to be defined on S. 
Standard electromagnetic theory can then be used to deter- 
mine, using the surface equivalence theorem, the compo- 
nents of the scattered electric field in the far field [27]. With 
these values, the distribution of the quantity 


E 


scatp 


Ie al + 


(0, b) = lim 4ar° 


ro 


a ei (19) 
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is computed, where (7, 0, @) denotes a standard spherical 
polar coordinate system. For the results presented here, the 
collection surface, S, is taken to be the surface formed in the 
mesh by removing all but the first layer of tetrahedra that are 
attached to the scattering surface and it is the distribution 
of the quantity 


RCS (0, 6) = 10 logio(y) (20) 


that is displayed. 


7 Basic Algorithm: Initial Validation 
7.1 Scattering by a2/, PEC Sphere 


The first example involves modelling the interaction between 
a monochromatic plane wave, of wavelength Ay, and a PEC 
sphere of diameter 2A). Simulations are performed initially 
using the standard FDTD algorithm and regular structured 
meshes with mesh spacing 6 = A)/15 and 6 = A,/60. 
Fig. 10 shows the staircased representation of the surface 
of the sphere on these structured meshes. The simula- 
tion is repeated using the unstructured FDTD algorithm, 
with a conforming unstructured mesh with global spacing 
6 = A) /15. For this example, cuts through a typical hybrid 
mesh are shown in Fig. 11. The analytical [26] and computed 
RCS distributions are compared in Fig. 12. It is apparent 
that the distributions computed on the fine structured mesh 
do not achieve the accuracy of the solution computed on 
the conformal unstructured mesh. This indicates that, for 
this problem, a significantly finer structured mesh must be 
used with the standard FDTD method to achieve accuracy 
comparable to that produced with the geometry conforming 
unstructured mesh FDTD method. 

The comparison is repeated for the case of 2A, dielectric 
sphere, with the relative values e, = 2 and y, = | assigned 
to the material parameters in the dielectric. In addition, it is 


(a) (b) 


Fig. 10 Representation of the surface of a 2A) PEC sphere using 
structured hexahedral meshes with spacing a 6 = A) /15; b 6 = A, /60 
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Fig.11 Cuts through a typical mesh employed for the analysis 
of scattering of a plane wave by a 2A, dielectric sphere: a view of 
the discretised dielectric sphere; b detail of a cut through the mesh 
showing the different cell types of tetrahedra (white and blue), pyra- 
mids (yellow) and hexahedra (red). (Color figure online) 


assumed that o = o,, = 0. The corresponding RCS distribu- 
tions are displayed in Fig. 13. In this case, the fine structured 
mesh is 64 times larger than the unstructured mesh and the 
corresponding allowable time step is four times smaller. 
The result is that, for this example, the fine structured mesh 
FDTD algorithm CPU requirement is 256 times larger than 
that of the unstructured FDTD method. 
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The same example is used to demonstrate the order 
of convergence of the method on a general unstructured 
mesh. Meshes with global spacing of Ap /4, Ap /8, Ap /15 
and A,/20 are employed and the L? norm of the errors in 
the computed distributions of the RCS are evaluated after 
20 cycles of the incident wave. Fig. 14a shows a conver- 
gence rate of 1.59 for the co-polarized RCS and Fig. 14b 
shows a corresponding rate of 1.48 for the cross-polarized 
RCS, giving an averaged value of 1.54 for the order of 
convergence. For a lossy dielectric sphere, with the same 
values for €,, uw, = land o,,, but with the electric conduc- 
tivity o = 0.7 S/m, the bi-static RCS distributions, com- 
puted after 20 cycles of the incident wave, are again seen 
to be in visually in a good agreement by comparing it with 
the analytical solution as can be seen in Fig. 15. 


7.2 Scattering by a Coated 2A, Sphere 


We now consider the scattering of a plane wave by a 
sphere, of electrical length A), that is coated with a uni- 
form dielectric, of thickness A). A cut through the mesh 
used to represent the sphere and the coating is shown in 
Fig. 16a. For the simulations reported here, the PML con- 
sists of 10 layers of hexahedral elements and the complete 
mesh consists of 876, 116 cells, with 1, 673, 527 Delaunay 
edges and 2, 076, 019 Voronoi edges. 
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(a) (b) 


Fig. 16 Scattering of a plane wave by a coated spherical object: a 
view of a cut through the mesh employed, showing the coating (blue 
cells) and the the object (red cells); b view of the electric scattered 
field E in the coating and on the object. (Color figure online) 


‘scat,y 


Initially, the inner sphere Fig. 16a is assigned the mate- 
rial parameter values ¢, = 1, uw, = 1, o = 10° S/m and 
o,, = 0 S/m, so that it effectively is a PEC. The coating is 
assumed to consist of a conducting dielectric material, with 
parameters €, = 2, uw, = 1, o =0.7 S/m and o,, = 0 S/m. 
The computed bistatic RCS distribution is compared with 
the analytical distribution in Fig. 17. An illustration of how 
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the scattered electric field propagates through the sphere is 
given in Fig. 16b. Finally, the sphere is considered to be a 
dielectric, characterized by the parameters e, = 2, w, = 1, 
o =0,, = 0S/m, and the dielectric coating is characterized 
by the parameters €, = 1.5, w, = 1, 0 =6,, = 0. The com- 
puted bi-static RCS distribution is compared with the ana- 
lytical distribution in Fig. 18. For these two examples, the 
computed and analytical RCS distributions are seen to be in 
excellent agreement. 


7.3 Scattering by a 15/, PEC Sphere 


To illustrate the performance of the procedure when it is 
applied to problems involving bodies of larger electrical 
length, we consider the simulation of the interaction between 
a plane single frequency incident wave and a PEC sphere of 
electrical length 15). The incident wave propagates in the x 
direction, the global element size is 6 = A)/10 and the outer 
surface of the PML region is, at its closest point, located at a 
distance of 2A, from the surface of the sphere. The staircase 
boundary is placed at a distance of A) /2 from the surface 
of the sphere and the thickness of the PML region is set 
equal to A,. The generated mesh has 7, 006, 329 vertices, of 
which 2, 460, 486 vertices are located within the PML. The 
unstructured mesh region near the sphere has 848, 205 ver- 
tices and, before the application of the mesh enhancement 
procedure, 4.3% of the cells in this region are such that the 
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centre of the cell circumsphere lies outside the cell. Follow- 
ing the application of mesh enhancement, the percentage of 
such elements is reduced to 1.1%. By allowing the merging 
of adjacent cells to form polyhedron cells, this percentage is 
further reduced to 0.2%. The majority of these undesirable 
cells are located adjacent to the staircase boundary. A view 
of the surface mesh and of a cut through the final mesh is 
shown in Fig. 19. 

The solution is advanced through 30 cycles of the inci- 
dent wave, using 180 time steps per cycle. A view of the dis- 
tribution of the computed contours of the E,,,,. component 
of the scattered electric field, on the scattering surface and 
on a cut through the mesh, is displayed in Fig. 20. For com- 
parison, the simulation is repeated using an explicit nodal 
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low order finite element time domain scheme (FETD) [28, 
29]. The FETD scheme is significantly less efficient requir- 
ing 8 times longer to advance the solution for one complete 
cycle of the incident wave. The RCS distributions computed 
with both schemes are seen to be in excellent agreement with 
the exact distribution in Fig. 21. 

To demonstrate the effectiveness of the advocated mesh 
generation approach, an alternative mesh is generated, using 
a basic Delaunay algorithm with automatic point inser- 
tion [13], in the region between the staircase and the sphere. 
In this case, 47% of the tetrahedra are such that the centre 
of the circumsphere of the tetrahedron is located outside 
the tetrahedron. This percentage is reduced to 4.3% after 
applying the traditional mesh enhancement procedures of 


Fig. 17 Scattering of a plane GoNolume 30 GoNolure = 
wave by a 2A, PEC sphere Analytic Analytic 
coated with a conducting _ a ] 
dielectric: a co-polarized RCS S S 
distribution; b cross-polarized 8 8 
RCS distribution 7a _ 
2 ae) 
g g 
Q QD 
a a 
-20 
0 20 40 60 80 100 120 140 160 180 0 20 40 60 80 100 120 140 160 180 
6 (deg) @ (deg) 
(a) (b) 
Fig. 18 Scattering of a plane 30 ease ee 30 Eau ae 
wave by a 2A) dielectric sphere Analytic Analytic 
coated with a dielectric: a co- 20 
polarized RCS distribution; b S S 
cross-polarized RCS distribu- 8 Q me 
tion x os ; 
g Z 9 
s s 
Fa & -10 
-20 
-20 1 1 1 1 1 1 1 1 1 1 i A A A 1 
0 20 40 60 80 100 120 140 160 180 0 20 40 60 80 100 120 140 160 180 
6 (deg) 6 (deg) 
(a) (b) 


Fig. 19 Scattering by a PEC 
sphere of electric length 15): 
a view of the discretisation of 
the scattering surface and the 
discretisation on a cut through 
the domain mesh; b detail of 
this view 
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Fig. 20 Scattering by a PEC sphere of electric length 1519: view 
of the computed distribution of the E,..,. component of the scat- 
tered electric field on the scattering surface and on a cut through the 
domain mesh 


edge swapping, edge collapse and Laplacian smoothing [15]. 
Using the FETD scheme on this mesh, it is found that 162 
time steps are required to advance the solution for one cycle. 
If the merging criterion is applied, stability of the FDTD 
scheme on the resulting mesh requires the use of 5491 time 
steps to advance the solution through one complete cycle 
of the incident wave. This shows quite clearly that meshing 
techniques that are only designed to produce a good quality 
Delaunay mesh will not, generally, be suitable for use with 
the unstructured mesh FDTD algorithm. 
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8 Basic Algorithm: Geometries of More 
Complex Shape 


To demonstrate its predictive capabilities, the described 
approach is now applied to problems involving more com- 
plex geometries. 


8.1 Scattering by a PEC Aircraft 


This example involves the interaction between a plane single 
frequency incident wave and a PEC full aircraft configura- 
tion. The wave frequency selected is such that the electrical 
length of the aircraft is 10 . 

The main axis of the aircraft is aligned with the x direc- 
tion and the incident wave propagates in this direction, with 
the wave impinging directly onto the nose of the aircraft. 
The staircase surface is located at a distance 1, /2 from the 
scatterer and the thickness of the PML is taken to be Ap. 
Two meshes, with an element size 6 = A)/15, are created. 
The first mesh is designed for use with the FETD method, 
with the region of tetrahedral elements generated using the 
standard Delaunay approach with automatic point inser- 
tion [13]. This mesh contains 4, 274, 438 vertices, of which 
2, 129, 120 vertices are located within the PML region. The 
second mesh, for use with the FDTD co-volume algorithm, 
is generated using the approach outlined above and contains 
4, 378, 740 vertices, with 2, 208, 512 vertices in the PML. 
For this second mesh, the percentage of tetrahedra having the 
property that the centre of a circumsphere of a tetrahedron is 
located outside the tetrahedron is reduced, from 16 to 2.2%, 
by using the proposed mesh enhancement procedure. This 
percentage is reduced further to 0.5% by using the advocated 
merging process. The majority of these undesirable elements 
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Fig. 21 Scattering by a PEC sphere of electric length 154): comparison between numerical and the exact RCS distributions: a co-polarized RCS 
distribution (plane @ = 0); b cross-polarized RCS distribution (plane @ = 2/2) 
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again lie in the vicinity of staircase interface. A view of the 
discretised aircraft surface is given in Fig. 22 and a detail of 
this view together with a view of the mesh on a cut through 
the domain is given in Fig. 22b. The solution is advanced 
for 40 cycles of the incident wave. The FETD computation 
on the first mesh takes 6.5 times longer per cycle compared 
to the FDTD scheme on the second mesh. It is interesting 
to note that the time required for the generation of the mesh 
for the co-volume scheme is 81 min. A detail of the com- 
puted distribution of contours of the E,.,,, component of 
the scattered electric field on the aircraft surface is given 
in Fig. 23. The RCS distributions computed by the FETD 
method and the FDTD algorithm on these meshes are com- 
pared in Fig. 24. Note that, in the context of scattering by a 
PEC dodecahedron using a completely ideal mesh, we have 
demonstrated previously [11] the convergence properties of 
the co-volume algorithm and its improved performance on 
a given mesh, for problems involving singularities, over that 
of the FETD method. This means that the discrepancies that 
are apparent here between the two RCS distributions for this 
example are probably due to the inability of the nodal finite 
element scheme to adequately represent the solution in the 
vicinity of geometrical singularities, without the use of addi- 
tional local mesh refinement [19]. 


8.2 Illumination of a Radome by a Narrow Band 
Pulse 


This example illustrates the calculation of the transmission 
efficiency of a dielectric radome. Radomes of this type are 
typically used to conceal aircraft communication or radar 
systems. The radome is modelled as half an ellipsoid, with 
a lateral radius of 0.5 m, a length of 1 m and a thickness of 
0.05 m. The ellipsoid is constructed of IparPlexiglass, with 
a relative electric permittivity e¢, = 2.59 and a conductivity 
o = 0.015 S/m. The 200 MHz pulse 
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Fig. 23 Scattering by a PEC aircraft configuration of electric length 
1049: view of the computed distribution of the E,.,,., component of 
the scattered electric field on the aircraft surface 


Egy t)= ek /oy /22° anton kek) (21) 


[Eq. (21)], linearly polarised in the y direction and propagat- 
ing towards the front of the radome in z direction, is used 
to illuminate the radome, where rt = 3.336 x 10~° s denotes 
the pulse width. An indication of the mesh employed can be 
obtained from Fig. 25a. The mesh has 20 degrees of freedom 
per wavelength and uses one layer of elements to represent 
the radome. The PML region is discretised using 10 layers 
of hexahedral elements and the smallest distance between 
the inner boundary of the PML and the radome surface cor- 
responds to 8 cells. The complete mesh consists of 228, 371 
cells, 537, 375 Delaunay edges and 587, 500 Voronoi edges. 
Fig. 25b shows computed contours of the scattered electric 
component E,... of the pulse on the portion of the mesh 


scat 
shown in Fig. 25a. To verify that one layer of elements is 
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Fig. 22 Scattering by a PEC aircraft configuration of electric length 104): a view of the discretised surface; b detail of the discretised surface 


and of a cut through the mesh 
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Fig. 24 Scattering by a PEC aircraft configuration of electric length 1049: comparison between the computed RCS distributions: a co-polarized 
RCS distribution (plane 6 = 0); b cross-polarized RCS distribution (plane @ = 2/2) 


(b) 


Fig. 25 Illumination of a radome showing: a a view of the coarse 
mesh used, with the yellow cells representing the dielectric and the 
white and blue cells representing free space; b a view of the com- 
puted component E,,,, of the scattered electric field. (Color figure 
online) 


sufficient to accurately model the dielectric, a new mesh is 
generated for the same radome, but using now 40 degrees of 
freedom per wavelength. The discretisation of the radome on 
the first mesh is shown in Fig. 26a and on the second mesh in 
Fig. 26b. The fine mesh contains 876, 314 cells, 1, 673, 283 


Fig. 26 Discretisation of the 
radome on a mesh using a 20 
degrees of freedom per wave- 
length; b 40 degrees of freedom 
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Delaunay edges and 2, 076, 114 Voronoi Edges. The fre- 
quency dependent transmission, 7T(f), which corresponds to 
a measure of the amplitude of the total electric field divided 
by the amplitude of the incident electric field, at a point rp 
behind the dielectric object is evaluated as 


FT [ErosQo; | 


T= 20lee,.||__ 
v 510 ep [Eme(o> | 


(22) 


where FT corresponds to the Fourier transform of the electric 
field. The value of the transmission obtained from the finer 
mesh is compared with that obtained from the coarser mesh 
in Fig. 27. The small differences between these results sug- 
gests that acceptable accuracy may be achieved here by using 
just one layer of cells to model the thickness of the radome. 

To study the effect of different frequencies on the same 
mesh, two narrowband pulses are considered. Both have a 
pulse width of 1 x 10~?s, with the first centred at 1 GHz 
and the second centred at 1.2 GHz. In the frequency range 
where the pulse spectra overlap, the transmission is expected 
to be the same. The predicted variation of the computed 
transmission efficiency with frequency, for both pulses, is 
plotted in Fig. 28. The mesh generation for this problem 
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Fig. 27 Comparison of the transmission of a 200 MHz pulse through 
a dielectric Radome, with pulse width 3.336 x 10~° s, evaluated on 
the coarse and fine meshes 


Transmission [dB] 


0.8 0.9 1 1.1 1.2 
Frequency [GHz] 


Fig. 28 Comparison of the computed transmission spectrum, through 
a dielectric radome, for two narrowband pulses, with one centred 
at 1 GHz and the other centered at 1.2 GHz, on the mesh with 40 
degrees of freedom per wavelength 


is challenging due to the rather thin thickness of the layer 
and extra care needs to be taken when collecting the total 
field history inside the radome. To reduce the influence of 
the dispersion error the incident field values in the update 
Eqs. (13) and (14) are not taken from the known function 
of Eq. (21) but are interpolated from a 1D FDTD code run- 
ning in parallel. This 1D code runs with the same time step 
and with a uniform mesh, with the mesh spacing equal to 
the spacing of the regular hexahedra in the actual problem. 


9 Modelling Anisotropic Lossy Materials 
For anisotropic materials, the electromagnetic parameters, 


such as permittivity, permeability and conductivity, vary in 
different directions, so that they must be treated as tensors. 
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Such materials offer many new and interesting possibili- 
ties in engineering, e.g. a thin anisotropic coating may sig- 
nificantly change the RCS distribution of an aircraft. With 
their additional mechanical strength and weight advantages, 
compared to metals, composite anisotropic materials, with 
applications initially limited to stealth aircrafts, satellites 
and space shuttles, are now part of everyday life [30]. In 
addition, anisotropy is at the heart of the new metamateri- 
als, which have electromagnetic properties that do not occur 
naturally, such as a negative index of refraction [31]. It is 
clearly important to extend the capabilities of the FDTD 
method to enable the modelling of problems involving such 
anisotropic materials. The approach that will be followed 
was initially introduced within the context of a total field for- 
mulation [30], but our unstructured mesh extension employs 
a scattered field formulation. 

For a three dimensional lossy anisotropic dielectric 
medium, the Laws of Ampére and Faraday are expressed, in 
a scattered field form, as 


Fi |< + 6" Dey hee Hoc dl 
0A 


A 


= J OE,,,, dA 
— | (€-€)/) Ot (23) 
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/ [s + Slt") Baca -dA = — ¢ ee ‘di 


0A 

(i 7) OH,,,, dA 
/ H— Ho Or 
A 


a ‘| On Hine “dA 


A 

(24) 
where D denotes the electric flux density and B the magnetic 
flux. In addition, 7 is the unit matrix, é is the electric per- 
mittivity tensor, # is the magnetic permeability tensor and 
6 and G,, are the electric and magnetic conductivity tensors 


respectively. The constitutive equations 
D=zE B= fH (25) 


define the relationships between the electric and magnetic 
fields and the electric flux density and the magnetic flux. 
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9.1 Unstructured Mesh Algorithm 


For the solution of this equation set, the unknown at the centre 
of the ith Delaunay ee illustrated in Fig. 3a, corresponds 
to the projection, (Dycar js E scar)» Of the scattered electric field 
onto the direction of the edge. The unknown at the centre of 
the jth Voronoi — illustrated in Fig. 3b, corresponds to the 
projection, (By arj, A scar)» of the scattered magnetic field onto 
the direction of that edge. Direct discretization of the Laws of 
Ampere and Faraday then leads to the equations 


Pav, By B'!/2 | 
2 seat jj 


_ Até, /T. Até 
eit fe 2] (fp 


J 
+ At — Pe - Fa, Hil; 


SCati;_ I,, 


j fl 
=(44= Ho ne -H’,|;] )-)) 
(26) 
and 
2 -1 = 
n+1 + ALiGnE x, 5 Ato yé 5, n 
Di = I + 2 P= 2 De icls 
MY 
Low ppti2w i, = peti? 
7 At AY » Fiche ©; — Og, Ei , 
spats tex 5, 0 n+1/2 ]).e) 
(Eqy — E01) Eine |} 08 
(27) 


Here, the bracket notation (F, é;) is used to denote the pro- 
jection of a field vector F along the ith edge. At interface 
boundaries, the material parameters in the Eqs. (27) and (26) 
are not constant and are replaced by averaged values é,,,, 
Agys Gq, and G,, . In this case, every component of the mate- 
rial parameter tensor is averaged following the approach 
employed earlier for the isotropic case [32]. To represent 
the vector-matrix multiplications that are required, we define 


. (2 Ree = fe Liew, 
gy = I + a) a, = i= 27 


(28) 
6 fae _ (Aten; 
Gy, =| I+ 5 a, =| 1- 5 

(29) 
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(31) 
and these definitions enable us to write Eqs. (26) and (27) as 


Bias = (ai (4y-Brea!"ly + Zal))-8)) (2) 
and 

Devan = (Aes (4e-Dycarli + Zoli) -€) (33) 
The quantities D” |, E E | a and B"!/ ‘|, H" |, represent 


scat inc 
y 
field vectors computed at iti centre of the ith Delaunay edge 


and the jth Voronoi edge respectively. 


9.2 Obtaining Approximated Field Vectors 
from Edge Projections 


To perform the necessary matrix—vector multiplications in 
Eqs. (32) and (33), we need to obtain the field vectors D,..,, 
and B,,.,, and associate them with Delaunay edge i and Voronoi 
edge j respectively. Unfortunately, we cannot get exact full 
field vector components directly from edge projections, but 
we are able to approximate the full field components at any 
location in the mesh. To achieve this, we assume that, with a 
set of three orthogonal vectors v,, V5, V3, a general vector x can 
be reconstructed as 


3 
=)'P, (34) 
i=1 


in terms of the projections 


aw y 24,953 

a 7 ra 2) 
As the mesh is unstructured, we cannot use this method 
directly to obtain an exact field vector (Please note that Ein- 
stein notations is not employed in this paper). The trivial 
solution would be to consider all the edges connected to 
one node, as depicted in Fig. 39b, and solve a system of 
equations. This will give us a field vector at node n,. As one 
edge is always formed by two points, we can do the same 
for point n, and average these field vectors and project to the 
corresponding edge. For example, to obtain the displacement 
field vector D,.,, on the edge i, denoted by the red line in 
Fig. 29a, we have to construct the set of equations 

PD... = (D é, (36) 


scat scat * e; ) 


where 
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Fig. 29 Interpolation of the field 
vector at node n,: a tetrahedra 
in the Delaunay mesh contain- 
ing the point n,; b edges in the 
Delaunay mesh containing the 
point n,. (Color figure online) 


a ej i, €;,%i, &; ei. 
P= Bei, i,@i, 7 ee (37) 
ee. e e;e 


and 6; is the normalized Delaunay edge vector correspond- 
ing to Delaunay edge i. The entries in the matrix P, based 
upon the three components of the vector €; and the projec- 
tions D,..,, : €;, are all known. The displacement field vector, 
Dy ai, at a node belonging to Delaunay edge i is approxi- 
mated by considering the sum of the system in Eq. (36) for 
each Delaunay edge connected to the node, as depicted in 
Fig. 29b. In this case, we solve the system 


N 
P’Dyoat = >, (Decang * ®q)&, (38) 


q=1 


N 3 3 
tim = Dy dy Dy Calan (39) 


Here, P’ and Dycatg * &, are known and N is the number of 
Delaunay edges connected to the node. This system of equa- 
tions is solved locally, node by node, resulting in an approxi- 
mated field vector at all nodes of the Delaunay mesh. The 
computed field vector is then linked to the corresponding 
edges. A Delaunay edge is assumed to connect the two nodes 
n, and n, and the displacement field vector, associated to the 
Delaunay edge i, is obtained by averaging the electric field 
vectors at nodes n, and 1, i.e. Deca.) = (Dicsa, AFD seaeg ) /2. 
The same procedure is applied for approximating the mag- 
netic flux vectors B,.,,, but using now the Voronoi edges. 
However, unlike the standard FDTD scheme, where we have 
one update equation for each component of the field vectors, 
the approximations of the full vector fields obtained using 
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this method are not good enough for time iterating the field 
components. It is found that error accumulation causes the 
algorithm to become unstable. We modified it in order to 
circumvent this difficulty. 


9.3 Local Orthogonal Unit Vectors 


It is the projection, Deca; = Dyca + €;, Of displacement field 
vector D,,.,, onto a Delaunay edge e,, that is updated in the 
solution process. For an isotropic material, the electric 
permittivity and the magnetic permeability are scalars, so 
that updating the fields only involves scalar multiplication 
between the field projections and the scalar material proper- 
ties. For an anisotropic material, matrix—vector multiplica- 
tions, between the material tensors (€, 4, 6,6,,,) and the fields 
(Dgcat> Bscat)» ave required. To perform these multiplications, 
we create two linearly independent vectors for each Delau- 
nay and Voronoi edge. Using the stabilized Gram—Schmidt 
orthonormalization procedure, we generate three orthonor- 
mal vectors, as illustrated in Fig. 30. The first vector e,, 
to which the two linearly independent vectors (e’,,e’;) are 
added, remains unchanged during the whole process. Each 
set of three orthogonal vectors represents one local coor- 
dinate system, leading to as many local systems as Delau- 
nay and Voronoi edges. We have already seen how we can 
reconstruct approximated field vector components for each 
local coordinate system, using the field projections of the 
surrounding Delaunay or Voronoi edges. The field vectors, 
obtained in this way, are projected onto the two orthogonal 
vectors forming the local frame. The first vector of each sub- 
set can be immediately updated using the projection equa- 
tion employed in the isotropic case. For this projection, no 
error is induced by the field averaging. 
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Fig. 30 Generating the orthonormal local coordinate system 


9.4 Coordinate Transformation 


The material tensors are expressed in the global reference 
frame formed by the orthonormal vectors X, ¥ and z. Three 
orthonormalized vectors %’, ¥’, 2’ form the basis of each 
local coordinate system. Here, the symbol * denotes a unit 
vector and the symbol ’ refers to vectors or vector compo- 
nents in the local coordinate system. The Jacobian matrix 
J, defined by 


Ox dy dz 

7 oy oy oy 

J=|— —— Ee 
Ox oy Oz ee) 
Ox oy Oz 


can be used to transform from the global (x, y, z) to a local 
frame. Each component of the Jacobian can be interpreted 
as an amplification factor, describing how one coordinate 
in a given reference frame stretches, shrinks or rotates with 
respect to another coordinate in another reference frame. In 
our case, the Jacobian is purely a rotation matrix Jp, which 
can be directly calculated as 


xX x9 x4 

= = “Af A Ay A A A 

Jr=|y-X y -y yz (41) 
a! -X a! 9 a! -% 


This rotation matrix has unit determinant, i.e. det J R=. 
Using a simplified set of Maxwells equations the form invar- 
iance [33, 34] is easiest explained in the following paragraph 
by expressing them in the local coordinate system, as 


[ede = g wat [won =f Bat 
Ot ot 


A’ oA’ A’ oA’ 
(42) 
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The electric and magnetic fields, in local and global coordi- 
nates, are related as 
E’(r') =) 'E() Hr’) =(0,) Hr) (43) 
while a general material parameter tensor, M , 1S transformed 
to the local coordinate system as 


ql 


M'(r’) = — 
det (Jp) 


(44) 


In this way, the coordinate transformation may be absorbed 
completely into the material properties. 


9.5 Time Updating 


The magnetic field is updated over the dual mesh at the half 
time step, using Eq. (33), and the electric field is updated 
over the primal mesh at the full time step, using Eq. (32). 
We will describe the updating of the electric field projection 
EY yy Phe magnetic field H”__, is updated similarly. The coor- 
dinate transformation applied to the terms a,,,a,_, defined 
in Eq. (28), and the term é, defined in Eq. (32), gives 


-1 


=) = = Ato,,€,, ST 

a, =Jp{ 1+ 5 ae (45) 
= 5 5 Atoe,, ay i 

a,_ =Jr\ I- 5) Ire (46) 
Eo = Ine Ig (47) 


These values are stored at the corresponding Delaunay edges 
before entering the time loop. Within the time iteration loop, 
Eq. (31), for each Delaunay edge é,, is calculated and stored. 
This is readily accomplished, as the magnetic components 
i Hto5 lv 


F : : M 
in the circulation term )} ‘ 


- : are available, from the 
Sl scat yix Jix 


previous iteration, and the complete vector EY is a known 
function. Before updating the projection D”_, the vectors 
DD” and Zp are calculated using Eq. (38). These vectors are 


scat 
projected in each of the three orthogonal directions of every 
local frame, e.g. A =é,, e. é, from Fig. 30, and the matrix 
vector multiplication is performed. In this way, Eq. (33) in 
the local frame takes the vector form 
p’”*! = I [a D” 


scat E+ E- scat 


+ Zi (48) 


The constitutive equation may then be used to obtain the 
electric field projection. In practice, we only have to consider 
the first line of é. a due to the fact that our data storage is 
based upon field projections along the edges and not field 
vectors. As a result, we can obtain 
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m+1 _ =-1 'n+1 =—1 'n+1 =—1 +1 
ee ~~ Ey scat; + Ei18™ sented + Eav.13” scat.e! (49) 
These values are then again used in to compute the curl in 
Eq. (30). 


10 Anisotropic Material Modelling 
Algorithm: Initial Validation 


A series of examples, involving scattering of an incident 
plane wave by an anisotropic sphere, is performed to vali- 
date the revised algorithm and to demonstrate its numerical 
performance. Validation of the implementation is achieved 
by comparing the computed RCS distributions with those 
obtained from the open source program Discrete Dipole 
Scattering (DDSCAT) [35], which is an implementation of 
the frequency domain discrete dipole approximation [36]. In 
each case, the incident wave, propagating in the x direction, 
has free space wavelength A, and the electrical length of the 
sphere is 2A). The mesh employed is illustrated in Fig. 31 
and has an average edge length of A,/20. The PML region 
is located at a minimum distance of A, from the scatterer and 
is discretised using 10 layers of hexahedra. The complete 
mesh consists of 876, 116 cells, 1, 673,527 Delaunay edges 
and 2, 076, 019 Voronoi edges [23]. It is worth noting that 
our unstructured mesh FDTD scheme only needs 8 points 
per wavelength to produce accurate results in free space [8]. 
However, inside a dielectric, the wavelength, Ap,,), is less 
than the free space wavelength, A). For the anisotropic case, 
it is estimated that a mesh spacing of A) /20 mesh should be 
sufficient here. Despite the fact that coarser meshes could 
have been used in some cases, we decided to use the same 
mesh for all the examples. 


(a) (b) 


Fig. 31 Scattering by a dielectric anisotropic sphere: a cut through 
the mesh used to represent the sphere; b contours of E,.,,, on a cut 
through the computational domain 
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10.1 Electrically Lossy Anisotropic Sphere 
For the first example, the material parameters 

13 0 0 0.3 0 0 
E=| 0 16 0 6=| 0 05 0 

0 O 2.0 0 O 0.7 

100 000 
H=|010 6, =|000 

001 000 (50) 


are employed and matrix multiplications are required in 
Eq. (27). Steady state conditions were attained after ten 
cycles of the incident wave. The computed RCS distribu- 
tions are shown to be in excellent agreement with the DDA 
results in Fig. 32. 


10.2 Magnetically Lossy Anisotropic Sphere 
In the next example, the material parameters of the sphere 


are characterised, by an anisotropic permeability tensor, 
together with magnetic conductivity, as 


100 000 
—é=1010 6=)}000 
001 000 
13 0 0 0.3 0 0 
H=} 0 160 6, =| 0 0.5 0 
0 0 2 0 O 0.7 (51) 


Now, matrix multiplication is required in Eq. (26) and 
steady state is again reached after ten cycles of the incident 
wave. The computed RCS distribution is compared with the 
DDA results in Fig. 33. 


10.3 Analysis of the Full Anisotropic Model 


To investigate the time penalty that results from the use of 
the anisotropic model, two simulations are performed on the 
same mesh. In the first example, the sphere is taken to be 
an isotropic lossy dielectric material. In the second exam- 
ple, the sphere is modelled as an anisotropic lossy dielectric 
material, characterised by the parameters 


13 0 0 0.3 0 0 
z=| 0 16 0 é=] 0 05 0 

0 0 2.0 0 0 07 

i2 0: 6 0.40 0 
a=| 0 14 0 &, =| 0 03 0 

0 0 18 0 0 05 (52) 
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Fig. 32 Scattering of a plane 
wave by a dielectric sphere of 
electrical length 2A, with aniso- 
tropic permittivity and electrical 
conductivity: a co-polarized 
RCS distribution; b cross-polar- 
ized RCS distribution 
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Fig. 33 Scattering of a plane 
wave by a dielectric sphere 
of electrical length 2A, with 
anisotropic permeability and 
magnetic conductivity: a co- 
polarized RCS distribution; b 
cross-polarized RCS distribu- 
tion 
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Comparison with DDA results is not possible for the second 
example, as the DDA code only allows for the modelling of 
non-magnetic materials. It is found that the computational 
cost for the simulation involving the anisotropic sphere is ten 
times the cost of the simulation for the isotropic sphere. This 
extra cost mainly arises from Eq. (38), which is a require- 
ment to solve a system of three equations for each Voronoi 
and each Delaunay node. The y and z components of the 
electric and magnetic fields respectively are represented in 
Fig. 34 

Given this additional cost, it is reasonable to ask if the 
current scheme is competitive when compared to the stand- 
ard FDTD method for this class of problems. We have previ- 
ously shown [23] that the accuracy of this scheme, with a 
Ao/15 unstructured mesh, compares favourable with that of 
the standard FDTD method, with a A) /90 or even a A) /120 
structured mesh for objects of arbitrary shape. This means 
that we get comparable results when using a mesh that is 
6-8 times coarser. A standard FDTD cell stores the elec- 
tric and magnetic field components at 12 edges and 6 faces 
respectively, corresponding to 18 degrees of freedom per 
wavelength. One FDTD cell may be discretised using 6 tetra- 
hedra, in the worst case scenario, leading to 19 edges and 14 
faces and a total of 33 degrees of freedom per wavelength. 
Although the number of degrees of freedom has nearly dou- 
bled, this is more than compensated by the fact that we may 
use a mesh that is, at least, 6 times coarser. Furthermore, for 
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the same volume in three dimensions, an unstructured mesh 
with 33 degrees of freedom produces the same accuracy as 
a structured mesh with 6° x 6 = 1296 degrees of freedom. 
The extra operations, needed for averaging and reconstruct- 
ing vectors, imply a cost penalty. For each node, we require 
18 averaging operations for each field, in addition to 3 x 3 
operations of vector reconstruction for each edge. For the 
6 tetrahedra lying inside one cube, the number of opera- 
tions will then be 2 x 8 x 18 +9 x 33 = 585. For one cell 
of a cartesian mesh, 8 operations are required for averaging 
the two offset components of a field vector, leading to 48 
extra operations on one cell node or 144 per FDTD cell. 
This implies that, when dealing with curved boundaries, 
to achieve the same level of accuracy, the standard FDTD 
scheme will actually require 6° x 48 = 10,368 operations 
per computational volume, whereas the current unstructured 
implementation has only 585 operations. 


11 Anisotropic Material Modelling 
Algorithms: Geometry of More Complex 
Shape 


To illustrate the predictive modelling capability of the 
enhanced solution algorithm, we consider again a problem 
involving the evaluation of the transmission efficiency of 
a pulse through a radome.The dimensions of the radome 
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(a) (b) 


Fig. 34 Scattering of a plane wave by a fully anisotropic dielectric 
sphere of electrical length 24): a contours of H,.,,, Shown on a cross- 
section through the mesh; b contours of £. 


scat,y SHOWN ON a Cross sec- 
tion through the mesh 


remain unaltered, but an anisotropic dielectric material, con- 
taining conducting fibres, is used to construct the radome. 
The mesh employed for the simulation, with 40 degrees of 
freedom per wavelength, is the mesh that was used earlier 
for the isotropic material case. The incident plane wave is 
again the narrow band pulse defined in Eq. (21), polarised 
in the y direction. 

Figure 35a is a schematic of a cut through a slab of 
composite material, in which the fibres are oriented in the 
y” direction. Reinforced carbon-carbon is such a type of 
composite where carbon fibres are reinforcing a graphite 
matrix. This material is for example used as heatshields 
for space shuttles. For the illustrated slab orientation, the 
material frame (X”’, 9,2’) and the global frame (X, ¥, Z) 
are identical. The situation is more complicated for the 
curved shell of the radome, as the orientation of the fibres 


Fibres Matrix 


(a) 
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changes in space. In this case, for each location on the shell, 
we have a material frame which, in general, differs from 
the global frame as illustrated schematically in Fig. 35b. 
Previously, we have only used the coordinate transforma- 
tion to pass from a global frame ((X, ¥, Z) coordinates) to a 
local frame ((X’, 9’, 2’) coordinates), linked to each Voro- 
noi or Delaunay edge. For a composite it becomes more 
complicated as we first have to make a coordinate trans- 
formation from the material frame, linked to the orienta- 
tion of the fibres (referred to as (X”, 9”, 2’’) coordinates) 
to the global frame ((X, Y, Z) coordinates). We then pass as 
before from the global frame, to a local frame (X’, y’, z’) 
which is linked to each Voronoi or Delaunay edge. Math- 
ematically these transformations are performed as follows: 
(8, 9,2") > Ip, > (RID) > Jpn > (8, 9,2’). Here, 
J rit = 1,2, are the two transformation matrices 


” 


: x! «x x” -y x" «7 
re, fat ” ”" " 
Sri =| y"-x y"-y y"-Z 
qx zl" -y Zw" 
7 x’ -x x’-y x’ -Z 
oe fee ! ! ! 
Jro =| y'-x yy yz 
z’-x Z-y ZZ 


(53) 


where Tay corresponds to the first transformation matrix 
describing the transformation from the material to the 
global frame, while Jp, corresponds to the second trans- 
formation matrix from the global to the local frame. This 
is therefore identical to the matrix Jp. Both Jp;,i = 1,2 are 
computed using Eq. (41), but with respect to the corre- 
sponding coordinate system leading to Eq. (53). We create 
an orthonormal material frame, which simplifies the coor- 
dinate transformation as the transformation matrix becomes 
a simple rotation matrix. The procedure here is illustrated 
in Fig. 36a. Initially, we only consider Voronoi edges at the 
dielectric interface. Each Voronoi edge (Vor1,blue) crossing 


(b) 


Fig. 35 Composite material: a orientation of fibres in the material frame; b orientation of the fibres inside the radome 
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Von 


(b) 


Fig. 36 Radome of composite material: a construction of the orthonormal material frame; b distance between cell centres and material frame 


located at the interface 


the face of a tetrahedron is surrounded by three Delaunay 
edges. We select two of these edges, Del, and Del,. By con- 
struction, these edges are parallel to the surface and taking 
the cross product between them the creation of the vector 
z"’ = Del, X Del, that is perpendicular to the surface. In the 
next step we assign x’ = Del,. By construction x” and z” 
are perpendicular to each other. To create the third vector 
y” perpendicular to x” and z”, the cross product between 
z” x Del, =z" Xx" =y” is taken. Following normaliza- 
tion, we have one orthonormal material frame (x”, y”’, 2”) 
for each face of the dielectric interface. In the next step, the 
material frame is linked to the cells inside the composite, by 
comparing the distances between the intersection points of 
a Voronoi edge with the interface and the circumcentre of 
all cells C; inside the dielectric. After finding the smallest 
distance, we link the local material frame (%’, 9’, 2’) from 
the surface to the corresponding cell. This process is illus- 
trated in Fig. 36b. To illustrate this process, consider a typi- 
cal material parameter tensor M, defined with respect to the 
material frame (X’’, 9”, z’’). The first coordinate transforma- 
tion converts the material parameters, from the material to 
the global frame, according to 


(54) 


The second coordinate transformation then converts the 
parameters from the global to the local frame, a process 
that has been described earlier. If we consider, for example, 
Eqs. (28) and (29), the material parameters €, 6, 4, G,,, have to 
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be replaced by €”,6”, j1”, 6" in each equation. The radome 
is characterized by the material parameters 


2.59 0 0 100 0 0 
é”=| 0 2.59 0 é” =| 0 100 0 

0 0 47 0 0 10 

100 000 
i’ =|010 &” =1000 

001 000 (55) 


Due to the orientation of the fibres, the properties 
of the composite are the same in the y” = mi” and the 
%” = directions, but differ along the 2” = fi” axis (see 
Fig. 36). This means that for the material parameters 
M(1, 1)" = M(2,2)” 4 M(3,3)", and the off-diagonal com- 
ponents are 0, for €”, «”’. The magnetic permeabiility tensor 
is set to the unit matrix I, i’ = 1 corresponding to free space 
and 6!” = 0. Typically, a material that minimally attenuates 
the electromagnetic signal transmitted or received by the 
antenna is used. 

The computed transmission for a 1 GHz narrowband pulse 
through the composite radome is represented in Fig. 37. 


12 Modelling of Bi-isotropic Materials 


Innovative fabrication techniques are currently being 
employed to create new chiral materials for use in communi- 
cation technologies. For these materials, a metallic structure, 
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Fig.37 Transmission of 1 GHz pulse through a radome constructed 
of a fiber based composite dielectric 


with the size of around 1/10 of the wavelength of interest, is 
embedded in a dielectric matrix. These inclusions need to 
be much smaller than the wavelength, so that the incoming 
wave sees the matrix and the inclusions as one homogeneous 
material. Often, these inclusions take the form of a combina- 
tion of a straight wire and a ring, to which the electric and 
magnetic field respectively will couple, leading to a material 
(Fig. 38a) with a frequency dependent response. For such 
materials, the electric permittivity and the magnetic perme- 
ability are frequency dependent parameters. Furthermore, 
for specific shapes or arrangement of the metallic inclu- 
sions, e.g. the use of left or right handed helices, as shown in 
Fig. 38b, the electric and magnetic fields will be related by 
an additional frequency dependent material parameter in the 
constitutive equation, referred to as the chirality. The chiral- 
ity allows for a negative index of refraction with, simultane- 
ously, a positive relative permittivity and permeabilty.Chiral 
materials can also create phenomena such as optical rotatory 
dispersion (ORD) and circular dichroism (CD). In ORD, 


Fig. 38 a Metamaterial devel- 
oped by Smith et al. [37], b 
chiral metamaterial consisting 
of left handed helices in an 
epoxy resin matrix 
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the polarization is continuously rotated inside the material, 
while CD refers to the change of polarisation of the incident 
wave, from linear to elliptical, due to the different absorption 
coefficients of a right and a left circularly polarized wave. 
Such materials are being used to create smaller antennas, 
super lenses [38], polarizers and radomes or RCS reducing 
materials. 

To extend the modelling approach to enable the simu- 
lation of problems involving chiral dispersive materials, 
a Lorentz model is employed for the electric permittivity 
and magnetic permeability and a Condon model is adopted 
for the chirality [39]. Several methods have been devel- 
oped for modelling frequency dependent materials within 
a time domain method, including the auxiliary differential 
equation method [40], the piecewise recursive convolution 
method [41] and the Z-transform technique [42]. As these 
three methods have been found to produce generally com- 
parable results, we decided to implement the Z-transform 
technique. This approach, which demonstrates good con- 
vergence close to the resonant frequency [43], is widely 
used in signal processing and can be viewed as the dis- 
crete version of the Laplace transform. It was first used to 
model chiral materials in 3D with FDTD [44], using the 
Z-transform to handle the frequency dependence of the 
material parameters. First order approximations were used 
in the Z-domain and the Z-transform coefficients were 
derived analytically. In a different approach [45], second 
order accuracy in the Z-domain is achieved by employing 
a bilinear transformation for calculating Padé approxim- 
ants [46]. In our approach, chiral materials are modelled 
by employing a generalization of this method. 


12.1 Problem Formulation 


It is not possible to employ a pure scattered field formu- 
lation when modelling problems involving frequency 
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dependent chiral materials. With such materials, an inter- 
mediate total field step is found to be required. In the 
frequency domain, the Laws of Ampére and Faraday are 
expressed in the integral form 


[Peale -dA = [Boe - dl 


AV oAV 


(56) 
+ [eo — Ey )jaE,,,.(@) -dA 
AY 
and 
JOB seat dA = — Ey cat - di 
AD oAP (57) 
+ [mo — Uy JOH jy. “dA 
AD 


where j = //—1. In the total field formulation, the corre- 
sponding equations are 


[Put -dA = J Ben . at | egioB (0 -dA 


AV oAV AY 

(58) 
and 
[ibe -dA = — ; pee -di+ / Mj OH ine “dA (59) 
AD oAP AP 


so that the terms €(@) and y(@) do not appear. Constitutive 
equations for chiral materials are expressed in the form 

jk JK 
Djo, = EE, + Mor Bio = MA: — Bor (60) 


where c is the speed of light and x the chiral parameter cou- 
pling the electric and magnetic fields. The Lorentzian model 


(€¢ ~ €c0) 


€(@) = €,, + ——_—__————__ 
() = Eg 2 +20, §,.jo — @ 
2 
(Me = Hoo )@ 
H(®) = Hoy + —" 
@, +20, 6,j)0-@ (61) 


is employed to determine the variation of the permittivity 
and the permeability with frequency. Here, the subscript co 
denotes high frequency values, while the subscript 7 refers 
to the low frequency limit. The values w,, @,, refer to the 
resonance frequency and €,, &,,, to the damping coefficients. 
The chiral parameter is obtained, from the Condon model, as 
2: 
jx) =@ = EE" 


a 62 
O, + 2a,& joo? (62) 
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where «, denotes the resonance frequency, while &, cor- 
responds to the damping coefficient for the Chiral material. 
The coupling constant, t,, defines the magnitude of the 
chirality. 

Using these relations, Eqs. (56) and (57) may then be 
expressed as 


[ie Em -dA = [Ben ‘di 


AV oAV 


+ | eiaB . aa- [J -dA (63) 


AY AY 
+ [Ja “dA 
AVY 
and 


[sete “dA = — / Eyeat -di+ Uj OH, “dA 


AD oAP AP 
- [In -dA — 38 -dA 
AD AD 
(64) 
where 
‘ joK 
Sin = jo(u _ Hoo Hot ich = ape H,,, (65) 
‘ jax 
Vee = jole ~ E oo) E ror Ie = _ E,,, (66) 


12.2 Discrete Equations and the Z-Transform 


To convert Eqs. (63) and (64) from the frequency to the time 
domain, we use the substitution ja > 2 and then apply the 
FDTD time stepping scheme to discretize the resulting equa- 
tions. The bilinear transformation 


oleae 
~ Atl+Z-! 


jo (67) 
is used, to convert a frequency dependent function to the 
Z-Domain [42]. Equations (65) and (66) can be expressed, 
with second order accuracy in the Z-Domain, by adopting the 
Padé approximants bi bi u ai, al for ij = ee, ke, hhor kh. 
In this manner, we obtain, for example, the approximation 


be + BEZ"! + BZ 


= (68) 
ee = = 
larZ 1 + ayZ 2 


Transfer to the time domain is direct, as Z~” acts as a 
backwards operator in time, e.g. for m = 1, Zz ye =I". 
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Rearranging Eq. (68) and, for an arbitrary field, F using the 
notation Z~”F”" = F”"”, we obtain 


Jit) = WwW" + bE! 
wrt bE! — 
ee 


ae’ n+1 
tot J 


+ BNE — 9S c. (69) 
Similar approximations may be written for the terms J, ” 
1 and i isla Using this notation, discretising the tile 
dependent form of Eqs. (63) and (64) on our unstructured 


mesh leads directly to the update equations 


1 
EB"! a [2&0 E° 
tot,t 2E oo + be At é tot,i 
1 MY nt+1/2 
i n+1/25V ine,i 7 
+ Ary 2) Ly Heats +e (70) 
ee 
and 
nt1/2 _ 1 n-1/2 
totj 2H a bh At tot.j 
MP nN 
—1 - n 'D inej 71 
+ Ats 2) AD y E scat. tigg + Ho at ( ) 


k=1 


ae Wing a i Sein \) 

The difficulty in these equations lies in the chiral material 
Cia and deg . For example, nei is 
required on the ae edges, bat it includes the electric 
field projections, that are only available on the Delaunay 
edges. The solution to this problem is illustrated in Fig. 39. 
For a tetrahedron, we have six Delaunay edges and the elec- 
tric field projection E,,i = 1,...,6 is stored on each edge. 
Using these projections, the electric field vector Ec, can be 
reconstructed and linked to the centre of the corresponding 
cell, in this case cell 1. To achieve this, the equation set 


coupling terms J,, 


Link reconstructed 
E vectors to cell center 


(a) 


Fig. 39 Determining the electric field on the Voronoi edges: a recon- 
struct the electric field vector from its projections on the Delaunay 
edges and link this vector to the corresponding cell center; b aver- 
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P Cell, ~ ell, 6; é; 
PE oon, = (Eceu, * &)€ (72) 


is constructed, where we again use the procedure described 
earlier. The resulting field vector is associated to the corre- 
sponding cell center. As a Voronoi edge connects two cells, 
we need to repeat this procedure for the second cell leading 
to E,,,,. The field vectors linked to these cells are then aver- 
aged and projected on to the corresponding Voronoi edge v 
according to (Eg, + Ecey)/2 - V, as illustrated in Fig. 39b. 
With the electric field defined on the Voronoi edges in this 
manner, we may then compute Wein: A similar process is 


adopted for the magnetic field, but in this case we have a 
point from the Delaunay mesh inside a Voronoi cell [32]. 


12.3 Solution Update Process 


To update the solution over one timestep, the total field for- 
mulation is employed as an intermediate step. In this man- 
ner, the solution is updated over one time step as follows: 


ye 


scat 


1. Magnetic field loop. Calculation of 


(a) reconstruct E,,,," from Delaunay edges and pro- 
ject to the Voronoi edges leading to E” 


scat, 
(b) add the incident field to obtain the total field vec- 
tor E” .= EE" + Et 
tot,j scaiy nc,j 
(c) compute J” att and pole 
: kt . mie n+1/2 
(d) compute HV) = Hy ae 


2. Electric field loop. Calculation of E”*!. 


scat, i 


(a) reconstruct Hf > from Voronoi edges and 


project to the Delaunay edges leading to oe 


scat,t 


(b) add the incident field to obtain the total field vec- 


tor yee hail aos 

ar n+1/24 oil ont 

(c) compute J Spa eee He and F . 
wit nel 

(d) compute fae i” ““tot,i inc,i 


Eceui 
Average and project 


v (Eceus - Eceu2) ;) 
2 


Eceu2 


(b) 


age the electric field vectors from the cell centers linked to the same 
Voronoi edge v and project them to the corresponding edge 
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12.4 Material Interfaces 


Boundary conditions at material interfaces are applied by using 
a weighted average process that is similar to that employed ear- 
lier for simulations involving isotropic and anisotropic media. 
However, for the Padé approximants, which are required in 
expressions such as Eq. (68), we only average the coefficients 
in the numerator. To illustrate this process, consider two tet- 
rahedral cells, separated by a boundary face, as illustrated in 
Fig. 40, and consider the calculation of the permeability for the 
corresponding Voronoi edge. Applying the bilinear transfor- 
mation to the frequency dependent component y,, — “4, from 
Eq. (61), leads to (u,, — P(z)), where 


P= 

" 1 +.a,z-1 + ayz~? (73) 
The magnetic flux is determined as 
B@); = H@H@); = He — P@))A@; (74) 


where the subscript i(= 1,2) refers to the medium 1 or 2. 
The average magnetic flux at the interface is calculated as 
L,B(z), + LBZ) 
BO)ay = > = 1 Hoo, — POYA®, 
ith 


+ W2(Hoo,2 — P(Z)2)H@)2 


where w,; = 1,/(, + /,) for i = 1,2 and J; corresponds to the 
length of the segment of the Voronoi edge between the inter- 
section with the Delaunay face and the center of the tetrahe- 
dron C; as shown in Fig. 40. 


Dielectric 1 


ae - / i ~ 
-- Voronoi Edge 
Delaunay Edges --7 “Length 1 — 
Intersection ~ _ 


at Boundary 


Dielectric 2 


Fig.40 Representation of a Voronoi edge and the corresponding cells 
at a material interface 
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13 Bi-isotropic Material Modelling: Initial 
Validation 


13.1 Rotation of the Plane of Polarization 


To verify the implementation of the scheme, the trans- 
mission coefficients for an electromagnetic wave incident 
upon a chiral slab are computed, along with an evaluation 
of the variation of the angle of the plane of polarisation 
of the incident field. The dielectric slab has dimensions 
0.3 m, 0.3 m and 0.1 m in the x, y, and z directions respec- 
tively and the mesh employed is such that 40 cells are used 
to discretise a length of 0.1 m. The electric permittivity 
and magnetic permeability of the slab are described by a 
Lorentz model and the chirality by a Condon model. The 
slab is located in free space, as depicted in Fig. 41. The 
slab has a real chirality, which induces an ORD of the 
incident field. In this case, the variation of the angle of 
polarisation, ®, can be computed as 


® = kool (76) 


where x is the chirality, w is the angular frequency of the 
incident wave and L is the thickness of the slab [47]. To 
compute this angle, we use a Fourier transform, which 
allows us to evaluate the variation of the amplitude without 
time dependence. A point in the Delaunay mesh is selected, 
such that, with respect to the incident wave, it lies behind the 
chiral slab. Then, for a wave propagating in z direction, the 
numerical value for the angle of polarisation is computed as 


Fig.41 Chiral slab (yellow) located in free space (blue). (Color figure 
online) 
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Fig.42 Rotation of the plane of polarization over time for 6 differ- 
ent values of the chirality. The black lines correspond to the theoreti- 
cal values and the coloured lines to the numerically computed values 
over several cycles. (Color figure online) 


FT (E,) 


Drum = ‘FT (E,) (77) 


Figure 42 shows the numerically computed rotation 
compared to the theoretical result, for different values 
of the chirality. For values of « less than 0.00407, the 
numerically computed angles of rotation of the plane of 
polarisation remain constant over time and are in good 
agreement with the theoretical values. For higher chirali- 
ties, however, the numerical solution shows very strong 
fluctuations indicating an instability of the scheme. We 
performed a stability analysis with respect to variations in 
the time step, the spatial step and the real and imaginary 
part of the chirality [48]. The results showed that the sta- 
bility is independent of the time step size. To investigate 
the effect of the spatial step size on the stability, a rotation 
per cell is defined by expressing the slab thickness, L, as 
L =nAz, where n is the number of cells through the slab, 
in the z direction, and Az is the edge length of a given cell. 
The cell normalised rotation of the plane of polarization 
is then defined as 
Py =< = Kevds (78) 
and the relative error in this quantity is computed. For 
a given real chirality, reducing the mesh size leads to a 
more stable scheme. A similar result is obtained if we use 
a given mesh and vary the chirality. Typically the lower 
the real part of the chirality the more stable the scheme. 
For a stable simulation and a real chirality, the results 
suggest that the scheme remains stable if ®,,,, < 0.4°. An 
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Fig. 43 for different damping factors €, = 0.01,0.025,0.5,0.1, the 
relative error of the computed angle over 20 cycles is plotted 


imaginary component of the chirality also has a stabiliz- 
ing effect on the scheme, as can be seen in Fig. 43. A 
coupling coefficient t, = 2 - 107!” s leads to strong insta- 
bilities in the scheme without damping. We increase the 
damping, using &, = 0.01, 0.025, 0.5, 0.1 in turn, and again 
plot the angles over 20 cycles. From the results, shown in 
Fig. 43, we see that the scheme is now stable well above 
the ®-,,, = 0.4° limit suggested for the non-damped case. 
It appears that, provided Im(x) > Re(x)/15, the scheme 
is stable, but, in this case, we are no longer dealing with 
a pure ORD. 


13.2 Reflection/Transmission on a Chiral Slab 


We consider the problem of simulating the interaction 
between an incident wave and a slab of bi-isotropic chi- 
ral material, with the object of determining the trans- 
mission and reflection coefficients. The chiral slab, 
which is assumed to be of infinite length and width, 
with a finite thickness L, fills the region 0 <z<L. The 
regions -co <z<0O and L<z<oware filled with free 
space. When the incident wave takes the form of a lin- 
early polarised narrowband pulse, as in Eq. (21), propa- 
gating in the z direction, the problem has an analytic 
solution [47]. 

In the numerical simulation, the thickness of 
the slab L=0.1 m is discretised with 40 cells and 
the chiral material is characterised by the val- 
ues €,=1.8€), €, =18€, w= 11g, Hy = 10M, 
O, = @, =@,=2HxX3.5 GHz, €=0.14, & =0.12, 
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é, =0.1 and t =1 ps. A point pl, in the free space 
region just ahead of the slab, is chosen for the compu- 
tation of the reflection coefficient, while the transmis- 
sion coefficient is calculated at a point p2, in the free 
space region slightly behind the slab. The co-polarized 
reflection coefficient, R,,, the co-polarized transmis- 
sion coefficient, T,,, and the cross polarized transmis- 


co? 


sion coefficient, T,,, are then determined as 


co? 


FT (Ey caty pl ) 
R.o recat a 
FT(E, 


inc,y,p1 ) 


= FT(E,-aty,p2) 
co T(E 


= F T(Escat.x,p2) 
“ FTE i,.5,2) 


(79) 


ineyy,p2) 


For a simple isotropic material, 7... is equal to zero. Good 
agreement is obtained in the comparison between the ana- 
lytical and computed distributions of the transmission coef- 
ficients shown in Figs. 44 and 45. The differences in the 
results arise because the analytical model assumes a slab 
which is infinite in both the x and y directions and, in addi- 
tion, the numerical dispersion inherent in the scattered field 
formulation leads to a shift in the peaks for the reflection 
coefficient. 


13.3 Chiral Sphere 


This example involves the simulation of the interaction 
between a plane wave and a chiral sphere. We consider a 
4/40 sphere and employ a mesh consisting of 4, 742, 854 
Delaunay cells, 1,727, 445 Delaunay vertices, 7, 604, 002 
Delaunay edges and 10, 544, 148 Voronoi edges, with 
10 layers of elements in the PML region. Three differ- 
ent values 0.8m, 1.0m and 1.2m are assigned to the free 
space wavelength A,. For Ay = 1.0m this corresponds to 
a sphere of electrical length 2A). The wavelength / inside 
a dielectric shrinks, with respect to free space, accord- 
ing to A=A,/Veé,u,—Kk?2. For the simulations, the 
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Fig.45 Comparison between analytical and numerical co-polarized 
reflection coefficients R.., 


values €, = 1.8€9, €,, = 1.6€, wy, = LL ao, Hy = 1.0 No, 
O, = @, = @, =2aX1GHz, € =0.12, € =0.12 are 
employed for the material parameters. A time step size 
At = 6.18 ps is used and the numerical solution is advanced 
through 20 cycles of the incident wave. The results produced 
for the RCS are compared to those produced in an analytical 
solution of this problem, using the program developed by 
Demir et al. [49]. 

The analytical and numerical RCS distributions are 
compared in Figs. 46, 47, 48, and 49 where, due to the 
symmetry of the RCS distributions, only the results for 
angles lying between 180° and 360° are shown. All the 
results demonstrate a very good agreement between the 
numerical and analytical solution for different values of the 
coupling constants 7, and the wavelength Ap. In particular, 
the very weak signal at high angles, from 300° to 360°, is 
well captured. 
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Fig. 46 Scattering of a plane 
wave by a chiral sphere 

with Ay = 1m, & = 0.3 and 

T, = 16.7 ps: a co-polarised 
RCS distribution; b cross-polar- 
ised RCS distribution 


Fig. 47 Scattering of a plane 
wave by a chiral sphere 

with Ay = 1m, & = 0.3 and 

T = 10 ps; a co-polarised RCS 
distribution; b cross-polarised 
RCS distribution 


Fig. 48 Scattering of a plane 
wave by a chiral sphere with 

Ay = 0.8m, &, = 0.3 and 

T = 16.7 ps: a co-polarised 
RCS distribution; b cross-polar- 
ised RCS distribution 


Fig. 49 Scattering of a plane 
wave by a chiral sphere with 

Ay = 1.2m, & = 0.3 and 

T, = 16. ps; a co-polarised RCS 
distribution; b cross-polarised 
RCS distribution 
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14 Conclusion 


We have demonstrated an unstructured mesh implementation 
of the classical Yee FDTD co-volume algorithm for compu- 
tational electromagnetics. The complexities associated with 
generating an appropriate primal Delaunay mesh, together 
with its Voronoi dual, have been identified and addressed. 
The approach was used initially for the simulation of wave 
scattering by PEC and isotropic dielectric objects, where 
the efficiency of the solution process, and the accuracy of 
the results, has been demonstrated. An advantage of the 
approach is that,with the same set of equations, it allows 
for the modelling of isotropic lossy dielectric materials 
or PECs by simply increasing the conductivities to very 
high values. This does not lead to any instabilities in the 
scheme.A description has been given of the extension of 
the approach, to enable the modelling of problems involv- 
ing anisotropic lossy materials and of problems involving 
frequency dependent bi-isotropic materials, with coupling 
between electric and magnetic fields. The incorporation of 
frequency dependent materials was accomplished by adopt- 
ing the Z-Transform method. 

This solution procedure must now be parallelised, so 
that its performance may be fully compared with that of 
other numerical techniques when applied to problems that 
are computationally more challenging. When this has been 
accomplished, the performance of the method can then be 
compared with that of the powerful current generation of 
hybrid finite difference time domain/finite element time 
domain methods [50-53]. 
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